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Abstract 

Plasmodium of true slime mold, Physarum polycephalum, is an amoeboid organ- 
ism, which spreads with developing tubular network structure and crawls on two- 
dimensional plane with oscillating the cell thickness. The plasmodium transforms its 
tubular network structure to adapt to the environment. To reveal the effect of the net- 
work structure on the oscillating behavior of the plasmodium, we constructed coupled 
map systems on two-dimensional weighted networks as models of the plasmodium, and 
investigated the relation between the distribution of weights on the network edges and 
the synchronization in the system. We found the probability that the system shows 
phase synchronization changes drastically with the weight distribution even if the total 
weight is constant. This implies the oscillating patterns observed in the plasmodium 
are controlled by the tube widths or cross-sections in the tubular networks. 

1 Introduction 

Plasmodium of true slime mold, Physarum polycephalum, is a unicellular amoeboid organ- 
ism. It spreads with developing tubular network structure and crawls on two-dimensional 
plane. Partial bodies of the plasmodium, whose cell thickness oscillate in period 1-2 min, 
interact with each other through protoplasmic streaming in the tubes [1]. The network 
structure of the tubes, which determines the interactions between oscillating partial bodies, 
has environmental dependency ||2j. For example, the plasmodium shows dendritic structure 
when the culture medium contains harmful chemicals, whereas it shows disk structure when 
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Figure 1: Morphology of plasmodium of true slime mold. Tubular network structures 
observed on agar medium with 1 mM KCl (left), and oat flakes (right). Bar, 1cm. 

the medium contains oat flakes, food for this organism (see Fig. [1]). From the viewpoint 
of the network theory [3l |3], these dendritic and disk structures can be characterized as 
tree-like and lattice-like networks, respectively [5]- 

Therefore, the plasmodium can be modeled as a two-dimensional network consisting of 
vertices that represent dynamical units, and edges that represent the interactions between 
the units [6]. In real plasmodium networks, widths of the tubes are not uniform but 
distributed largely. Because the width of the tube can be considered as the coupling 
strength in coupled biological oscillators of Physarum plasmodium [Tj [8], each edge in 
a Plasmodium network has its own "weight" that is directly connected to the coupling 
strength in the corresponding coupled dynamical systems. Thus, dynamical oscillation 
patterns observed in Physarum are thought to be controlled not only by the structure of 
the network, but also by the weight distribution on a network. 

Recently we investigated the latter effect by varying the weighting rules on a fixed net- 
work topology; this topology was designed as an intermixture of tree and lattice so that we 
can express observed plasmodial tube networks [6j. On these generated two-dimensional 
weighted networks, coupled phase oscillators are constructed and the propensity for collec- 
tive behavior, such as phase-locking and synchronous oscillation patterns, was investigated. 
As a result, we found that the distribution of edge weights in the networks strongly affects 
the global synchronization and spatiotemporal oscillation patterns, even if the network 
topology and the total weight are fixed. 

Can we apply this result to other coupled dynamical systems in general? As a step 
toward answering this question, in the present paper, we construct coupled logistic maps. 
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Figure 2: Designed topology with b = 2 and dmax = 3. 



instead of phase oscillators, on two-dimensional weighted networks, and investigate their 
dynamics in detail. 

2 Model 

2.1 Weighted planer networks 

To express the environment-dependent plasmodial tube networks, two-dimensional weighted 
networks whose topology is an intermixture of tree and lattice are constructed as described 
in ref.[6]. 

A network topology that includes a tree as well as a lattice is generated as follows. At 
first, a Cayley-tree with the branching number of b and the maximum depth of dmax is 
generated. The depth, d, of a given vertex is defined as the shortest path length from the 
tree's root to the vertex. Thus the depth of the root vertex is zero {d = 0). The number 
of vertices in depth d is given by n{d) = (6 -|- 1) • b'^~^ when d > 1. Thus this tree has 
N = 1 + Yld=T "(^) = ! + (& + - l)/{b - 1) vertices, and iV - 1 edges. Note that 

the degree (connectivity) of non-boundary vertex (i.e. vertices in the depth d / dmax) is 
6 -|- 1 in this tree. Next, the vertices in the same depth are connected side by side to make 
a ring. In this process, — 1 edges are newly added, and all vertices but the root increase 
their degrees by two. The resulting topology with b = 2 and dmax = 3 is shown in Fig. [2l 

Each edge on the designed topology is weighted using simple symmetrical rules with 
a small number of parameters. Firstly, we separate the edges into two groups. The first 
group consists of A^ — 1 edges originally belonging to the Cayley-tree, and the second group 
consists of the remainder A — 1 edges newly added to the tree. Below, weights on the 
edges belonging to the first and the second groups are denoted by w and v, respectively. 
Secondly, we assume the weights in a same depth are identical. Accordingly, the weights 
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can be determined by a single parameter d, such that w{d) and v{d). We define w{d) as the 
weight on the edge that connects a vertex in depth d—1 and one in depth d. Similarly v{d) 
is defined as the weight on the edge that connects vertices in the same depth d. Thirdly, we 
define the following parameters; R = w{d+ l)/'w{d), S = v{d)/w{d), and the total weights 
wtot given by 

^«tot= 5^n(d)K(i)+t;(d)]. (1) 
d=l 

For given R,S, and wtot, all of the weights on the edges in the network with the 
branching number b and the Illcixillllini depth, drnax 

are uniquely determined as 



w{d) 



+ + - 1]' 

WtotR' 

dmax{l+b){l + S)' 



d-i- ' ■ (2) 

hR = 1, 



and v{d) = Sw^d), for d = 1, 2, . . . , dmax- In Fig. [3l we illustrate how the ratios R and 5 
affect the appearance of one of the designed networks when wtot is fixed. The parameter 
R determines whether the network is center weighted or periphery weighted; the networks 
becomes center (periphery) weighted ones when R < 1 (R > 1). Another parameter S 
determines whether the network is tree-like or ring-like; the networks becomes tree-like 
(ring-like) ones when 5" < 1 (5 > 1). 



2.2 Coupled maps 

On the weighted, networks characterized by b^dmax-! 

R,S, and wtot, coupled maps are con- 
structed by putting a map on each vertex of the network and coupling the nearest neighbors. 
Dynamics of the coupled maps are given by 

Xi{t + 1) = fi[xi{t)] + ^ Ai,{/,[x,(t)] - fi[xi{t)]}, (3) 

j 

where Xi{t) is the state of the vertex i(= 1,2,... , N) at discrete time t, the function fi 
expresses a nonlinear map that generally depends on the index i of the vertex, and the 
matrix component Aij expresses the coupling strength between the vertices i and j. We 
assume the coupling strength is directly proportional to the weights on the corresponding 
edge. Here we set the proportional constant to be 1. Thus, if vertices i and j are connected, 
Aij{= Aji) is given by (i) Aij = w{d), when the corresponding edge e is in the first group 
of edges and the depth of the deeper vertex {i or j) is d, or (ii) Aij = v{d)^ when e is in 
the second group and the depth of the both vertices is d. 
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Figure 3: Effect of the ratios R and S on the appearance of the weighted network, b 
and dmax = 3 are used. 



In this paper, we used the logistic maps as dynamical units of the system; 

fi{x) = aix{l - x), (4) 

where Oj is a dynamical parameter that governs the uncoupled map dynamics. Since the 
amoeboid organism Physarum shows basically periodic behavior, we mainly consider the 
case when all of the dynamical parameters ai{i = 1,2, . . . , N) of the logistic maps are 
2 + Vq/2 ~ 3.22, the center of the region [3, 1 + \/6], where uncoupled maps show stable 
period- 2 solution. 



2.3 Phase synchronization 

To estimate the propensity for collective behavior, we calculate the probability that the sys- 
tem shows phase synchronization {Pgync), which is defined and calculated as follows. After 
weighting on each edge in the designed topology, we start calculation of Eq.(l3]) with initial 
conditions, Xi{0),i = 1,2, . . . , N, given randomly from [0,1]. After tj = 3000 iterations, 
we obtain the next Tq = 1000 values for each i; i.e. Xj(ry + l),Xi(r/ + 2), . . . ,Xi{Tf + Tq) 
for i = 1,2, . . . , N . Next, we calculate the phase distance between all pairs of vertices [5]. 
The phase distance dij between vertices i and j is defined as 



where is the number of times the observed state values of vertex i, Xi{t), show local 
minima during the time interval of [rj + 1,tj + Tq], and ^ij is the number of times the 
minima of Xi{t) and Xj{t) match with each other. When dij = 0, we say vertices i and j 
are phase synchronized [10]. If dij = for all pairs of vertices in a cluster of vertices, this 
cluster is called a phase synchronized cluster. If the size of such a cluster is identical to N , 
we say the whole system is phase synchronized. 

For each initial condition, we obtain the size of the largest cluster existed in the system 
and calculate the probability that the size is equal to the size of the system N . This prob- 
ability is identical with the probability Psync that the system shows phase synchronization. 

Because the coupling term in Eq.Q does not normalized with the strength = 
"^2^=1 ^ij, the logistic map may take an invalid value that is outside of [0, 1]. In this case, 
we stopped the calculation and started a new calculation with a new initial condition. We 

count these invalid cases as data showing no phase synchronization for the calculation of 
p 

^ sync- 
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Figure 4: (a) Probability that the coupled maps on a weighted tree shows phase synchro- 
nization {Psync) calculated in the R-W space. Each probability for a given set of R and 
W is given by the number of samples showing the phase synchronization divided by the 
number of observations (n = 100). (b) Pgync plotted against R for the cases of W =0.04 
(open circles), 0.21 (open squares), 0.34 (closed diamonds), and 0.45 (closed triangles). 



3 Results 

3.1 Coupled maps on weighted trees 

At first, we restrict our attention to the weighted trees, the case when there is no connection 
between vertices in the same depth, or 5 = 0. Here we consider the coupled identical maps 
with ai = 2 + y/Q/2 (for all i) on a Cayley tree with 6 = 2 and dmax = 3. Since S is fixed 
to be zero, only two parameters R and W determine the weight distribution on the tree. 
The mean weight W is given by = wtot/{^ ~ 1) = "Wtot/^l, because the tree has N = 22 
vertices and — 1 = 21 edges. 

We calculate the probability Pgync that the system shows phase synchronization in 
the space of R and W [Fig. [Ula)]. This figure shows the dependence of Psync on weight 
distribution (i?) strongly depends on the mean weight (1^). In cases of small mean weights 
[e.g. W = 0.04 shown in Fig. [Hb)], the coupled maps rarely synchronize because the 
coupling strengths are too small. In cases of moderate mean weights, Psync has a value 
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larger than 0.5 in the appropriate region of R. When W is relatively small, e.g. W = 0.21 
shown in Fig. |4)[|b), Psync has its largest value at the appropriate R less than 1. This implies 
that the coupled map system on a center- weighted network (R < 1) can phase synchronize 
easier than uniformly weighted {R = 1) or out-of-center weighted {R > 1) trees. On the 
other hand, when W is relatively large, e.g. W = 0.34 shown in Fig. \Mi>), Psync has 
its largest value at the appropriate R between 1 and 1.5. This implies that uniformly or 
slightly out-of-center weighted trees are appropriate for the synchronization of the coupled 
maps. In cases of large mean weights, e.g. when W = 0.45 shown in Fig. Hl^b), Pgync has 
very small value regardless of R because the maps frequently take invalid values outside of 
[0, 1]. With the same reason, Pgync is almost zero in the region near the lower left corner 
of Fig. IHa), where W is large and R is small. 

We also calculate Pgync in cases of randomly weighted trees, where each edge on a tree 
is weighted with a value given randomly from [W — A,W + A] (data not shown). As a 
result, we found that the more randomly weighted on the edges of a tree, the less probable 
the coupled map system on the weighted tree is phase synchronized. This implies that the 
case of A = 0, corresponding to a uniformly weighted tree, has Psync larger than any other 
randomly weighted cases at a given W. Note, however, that even in the case of A = 0, 
Psync is smaller than the largest Psync obtained in the cases weighted orderly with the 
appropriate parameter R. 

3.2 Coupled maps on the designed networks 

In the cases when there are connections between vertices in the same depth, or the cases 
of S* 7^ 0, we consider the coupled identical maps with ai = 2 + (for all i) on the 

designed network with 6 = 2 and dmax = 3. The probability that the system shows phase 
synchronization, Psync, is calculated in the parameter space of R and S for a given total 
weight wtot (see Fig. [5]). When wtot is small, for example when wtot = 1-0, the system rarely 
phase synchronizes regardless of the parameters [see Fig. [5Ka)]. On the other hand, when 
Wtot is large, for example when wtot = 4.5 or 10, the system always phase synchronizes 
regardless of the parameters, as long as the logistic maps take valid values, i.e. [0, 1] 
[see Figs. and (d)]. When the total weight is moderate, for example, in the case of 
Wtot = 2.0, Psync depends on the weight distribution determined by R and S [see Fig. [5{b)]. 
As shown in Fig. [5)[|b), Psync becomes very large when S > 0.5 and 0.6 < R < 1.0, implying 
that the coupled maps on center- weighted {R < 1.0) lattice-like {S / 0) networks are easy 
to phase-synchronize. 

If the dynamical parameter aj of maps are not identical but varied uniformly between 
1.8-F\/6/2 ~ 3.02 and 2.2-h\/6/2 ~ 3.42, the effect of weight distribution on the propensity 
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Figure 5: Probability that the coupled maps on the two-dimensional weighted designed 
network show the phase synchronization (Pgync) calculated in the R-S space when the 
identical dynamical parameter = 2 + (for all i) is used for each logistic map. The 

total weight (wtot) used for the calculation is shown in lower right corner in each figure. 
The number of observation for each parameter set of R and S is 100. 6 = 2 and dmax = 3 
are used. 
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Figure 6: [(a) and (b)] Probability that the coupled maps on the two-dimensional weighted 
designed network (Psync) calculated in the R-S space when the dynamical parameter 
for each map is given randomly from the region [1.8 + V^/2,2.2 + V6/2] . The total 
weight {wtot) used for the calculation is shown in lower right corner in each figure. Other 
parameters and calculation conditions are the same shown in Fig. [5l 

for phase synchronization changes. Figures Ella) and (b) show the probabilities that such 
coupled maps show phase synchronization (Psync) in the R-S plane when the total weight 
is 2.0 and 4.5, respectively. In the case of wtot = 2.0, the probability obtained in the system 
composed of non-identical maps has its peak at ~ 1 and 5 ~ 1, corresponding to the 
uniformly weighted network, although the probability is smaller than that obtained in the 
system composed of identical maps [Fig E^b)] for any given parameter set of R and S. 

4 Sufficient conditions that the system takes valid values 

As mentioned above, when wtot is large enough, the system always phase synchronizes as 
long as the maps take valid values, i.e. [0, 1] for the logistic map. In those cases, the region 
of the weight distribution {R,S), where the coupled map systems defined by Eq.([3|) are 
phase synchronized, can be derived by the sufficient conditions for the maps to take valid 
values. Such conditions are expressed mathematically as max[xi{t + 1)] < 1 (for all i), 
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where max[-] represents the maximum value. 
Since Eq.(l3]) can be rewritten as 



x,{t + 1) = (1 - s,)n + J2 ^ii/i' (6) 

j 

where Si = Aij is the strength of vertex i, we have 

max + 1) = < . (7) 

i(l-s,)/r" + E,^../r"' **>1' 

where we define the range of the function fi as [/™", f^"-^]. When the ranges of all of the 
maps are identical to ^e have 

^ ' . ' (8) 

+ >1. 

In the cases using the logistic maps, we have = and f™-"'^ = /(1/2) = a/4, where 
a is the dynamical parameter. While max[xi(i + 1)] = /"^"^ = a/4 is less than unity for 
any given a(< 4) when Sj < 1, max[xj(i + 1)] will be less than unity if 

imin ^ 

Si < = - (9) 

— jmax Jmin ^ ^ ' 

is satisfied when Si > 1. Therefore, the sufficient conditions that each map in the coupled 
identical logistic maps with dynamical parameter a takes values within the domain [0, 1] 
are 

-asi < 1, i = 1,2,... ,iV. (10) 

If we use the weighting rules mentioned above, the strengths are uniquely determined 
by the depth of the vertex, d. Thus, the strength can be written as a function of d, i.e. 
s{d). This function is given by 



sid) 



{b+l)w{l), d = 0, 

{l + 2S + bR)w{d), d=l,2,...,dmax-l, (11) 
(1 + 2S)w[djYiax) 1 d = draax- 



In those equation, w{d)^d = 1,2, ... ,dmax are given as Eq.Q. Using this function, the 
sufficient conditions can be rewritten as 

^as(d) < 1, d = 0,l,...,dmax (12) 
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Figure 7: (a) The region in which all of the logistic maps on the weighted tree never take 
invalid values when 6 = 2 and dmax = 3, 5, 10, and 30 [Eqs. ([13]) and ^^]. Each region 
includes the upper-right corner, (b) The same region when dmax = 3 and b =2, 3, 4, 5, and 
10. Each region includes the upper-right corner, (c) The region in which all of the maps 
on the weighted designed network never take invalid values when 6 = 2 and dmax = 3 [Eqs. 
(jl8p - (j2ip ]. Each region includes the lower-right corner. Numbers on/near the solid lines 
represent {a/4:)wtot- 



4.1 The case of weighted trees 

In the case of the weighted tree (i.e. the case of 5 = 0), the largest strength in the vertices 
is s(0) when R < 1, and s{dmax ~ 1) when i? > 1. Thus the conditions p2|) are respectively 
given by as(0)/4 < 1 and as{dmax — l)/4 < 1- Substituting Eq. ifTT]) into these equations, 
we have the sufficient conditions as follows. 

4(6- 1) [(6/2)'^™- - 11 , , , 

S a(6;i)(M-l)y.-.-l) - "-^^ (^^) 

< 4(>- l)l(tfi)-^" -l| „,^„ „ ^ J (1,) 

In the special case of bR = 1, where < 1 is satisfied because 6 is an integer larger than 
1, the condition is given by 

W< "^dmaxib-l) 

~ a(6-h 1)(6^™- - 1) ^ ^ 

Figure UHa) shows the region in which all of the maps never take invalid values when 
6 = 2 and dmax = 3, 5, 10, and 30. The regions include the upper-right corner {R = 3, 
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W = 0) of the figure. As the maximum depth increases, the region shrinks to that in the 
hmiting condition of dmax — > oo; 

W = 0, when R < 1 (16) 

Figure El^b) shows the same region when dmax = 3 and b =2, 3, 4, 5, and 10. The 
regions include the upper-right corner (R = 3, W = 0) of the figure. As the branching 
number increases, the region also shrinks to the upper. 

4.2 The case of weighted designed networks 

In the case of 5 7^ 0, the sufficient conditions that each map in the coupled logistic maps 
takes values within the domain [0, 1] [i.e. Eq. (jlOI) ] are given by 



1 
1 

4' 



< 


1, 


when 


< 


1, 


when 


< 


1, 


when 


< 


1, 


when 



(18) 
(19) 



^as{dmax-l) < 1, when 1< R < {I + 2S)/{1 + 2S - b) (20) 
1 . , . . , ^ 1 + 25 



(> 1) (21) 



where s{d),d = 0,1, dmax — 1) and dmax are given by Eq. (fTT|) . Figure El^c) shows the 
boundaries of the region in which all of the maps never take invalid values when b = 2 and 
dmax = 3. The region includes the lower-right corner {R = S = 2). Numbers on or near 
the boundaries represent {a/4:)wtot- These boundaries can predict the region in which the 
system shows phase synchronization with high probability when wtot is large enough. For 
example, the region of Pgync = 1-0 in the case of wtot = 10 [Fig. [5]^D)] is similar to the 
region in which Eq.([l2D is satisfied when {a/A)wtot = [(2 + \/6/2)/4] • 10 ~ 8 [see Fig.E^c)]. 

In the case of wtot = 2.0 and a = 2 + VQ/2, used in the calculation shown in Fig. EJb), 
the line of {a/4:)wtot ~ 1-6 corresponds to the boundary. Thus the all region without the 
small region including the upper-left corner {R = S = 0) satisfies the sufficient condition 
given by Eq. (fT^ [see Fig. EU^c)]. Therefore the dependency of the weight distribution on 
the probability Psync shown in Fig. EJ^b) does not include the effect that the maps take 
invalid values. 
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5 Discussion 



For many years, synchronization in coupled dynamical systems on complex networks has 
been studied by many researchers (see ref.|llj and references therein). Many of the studies 
considered a system composed of identical dynamical units on unweighted networks [12\ 
I14j . To elucidate synchronous phenomena observed in real networks [15], studies on 
synchronization of the systems on weighted networks are needed. Previous studies on 
coupled dynamical systems on weighted networks are based on the assumption that vertices 
in the networks have large degrees [16j . However, degrees of the vertices in many real 
networks, including tubular networks developed by the Plasmodium (Fig. [1]), are not so 
large. Therefore, in the present study, we have constructed coupled logistic maps on two- 
dimensional weighted networks whose vertices have small degrees (< b + 3). Using this 
system, we have shown that the probability that the system shows phase synchronization 
depends strongly on the weight distribution when the total weight is moderate [see Fig. 

We have also studied dynamics of phase oscillators on two-dimensional weighted net- 
works [6j, whose topology and weighting rules are the same shown in the present study. 
The probability that the system shows phase locking was calculated in the parameter space 
of R and S using a moderate total weight wtot- As a result, it is found that the high proba- 
bility region of phase locking (Fig. 4 in ref. [6] ) resembles the one of phase synchronization 
when the dynamical parameter for each map is randomly given (Fig. [6]). Therefore 
we speculate that the relation between probability showing collective behavior and local 
coupling strength, or weight distribution in a given network, is less dependent on types of 
dynamical systems, i.e., whether the system is composed of phase oscillators or maps. 

In conclusion, we studied the dynamics of coupled logistic maps on two-dimensional 
weighted networks as models for plasmodium of true slime mold, Physarum polycephalum, 
spreading with developing tubular networks. We found the weight distribution among edges 
of the networks strongly affects whether the coupled maps are phase synchronized even if 
the total weight is constant. This implies that synchronization in Physarum networks is 
controlled by changes in the width, or the cross sectional area, of each plasmodial tube 
that determines the coupling strength between partial bodies of the system. 
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